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Abstract. In antecedent experiments the electron screening energies of the d+d reactions in metallic en- 
vironments have been determined to be enhanced by an order of magnitude in comparison to the case of 
gaseous deuterium targets. The analytical models describing averaged material properties have not been 
able to explain the experimental results so far. Therefore, a first effort has been undertaken to simulate 
the dynamics of reacting deuterons in a metallic lattice by means of an ab-initio Hartree-Fock calculation 
of the total electrostatic force between the lattice and the successively approaching deuterons via path 
integration. The calculations have been performed for Li and Ta, clearly showing a migration of electrons 
from host metallic to the deuterium atoms. However, in order to avoid more of the necessary simplifications 
in the model the utilization of a massive parallel supercomputer would be required. 
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theories and models of atomic and molecular collisions and interactions - 26.20.+f Hydrostatic stellar 
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1 Introduction 

Following our first observation [TJHIIS] of the grossly en- 
hanced screening effect for the d+d reactions in metals 
resulting in an exponential-like increase of the astrophys- 
ical S-factors towards lower energies other groups repro- 
duced our results [4,5,6,7,8,9]. The measured screening 
energies are one order of magnitude larger than for D 2 -gas 
targets [10]. However, particular care is required for the in- 
terpretation of the experimental data because the special 
physico-chemical properties of the hydrogen compounds 
and the beam induced chemical reactions at the target 
heavily influence the obtained screening energy results O 
ITT] . Originally, the screening energy was extracted using 
the simple model of jT2] treating it as an shift of the kinetic 
energy. Thus, the screening energy is merely a parameter 
which describes the modification of the Coulomb barrier 
and not a real energy gain for the projectile. In a micro- 
scopic view it is universally valid that the screening effect 
depends on the impact of the electronic configuration of 
the environment on the Coulomb barrier of the entrance 
channel only, i.e. the pure Coulomb energy is simply mod- 
ified by a Yukawa factor W(r) = _J_ z p^* e e ~*Awith Aa 
being the screening length. As such the inferred screening 
energy is merely the second term in a Taylor-expansion 



of W(r), i.e. U e = j^- p x ^ e , and a coarse mathematical 
parametrization in the simple model. 

The experimental screening energy values could not 
be described by theoretical approaches to this extent. In 
[13]114] we presented an analytical model based on the di- 
electric function theory for the description of the screening 
effect. Its results, however, are below the measured screen- 
ing energies by a factor of 2 though better than other ap- 
proaches cited therein, e.g. [T5pl . Therefore, a first effort 
was undertaken to simulate the pre-reaction impact with 
an ab-initio quantum mechanical Hartree-Fock calculation 
which is able to consider the actual crystal structure while 
the analytical model only operates on averaged material 
properties. The concrete form of the electron density dis- 
tribution around the nucleus indeed influences the screen- 
ing energy to a great extent as was shown in [17] on a D2 
molecule with a time dependent Hartree-Fock calculation. 
Due to the many electrons involved in the simulation of 
the crystal the time dependency can not be retained on 
a workstation class computer. The model is described in 
the next section permitting the assessment of the neces- 
sary simplifications. 
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2 Mathematical Model 

The mathematical model, the simulation is based on, re- 
sults from a hierarchy of simplifications in the theoretical 
description of the physical reality leading to more or less 
severe deviations in the obtained values of the observables. 

The originator for it is the non-relativistic spin inde- 
pendent Schrodinger equation 
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containing both kinetic and potential operators for the K 
nuclei (0) and the N electrons (0 [jS Therein Z a is the 
charge and M a the mass of the nucleus a. The distances 
in the denominators of the three potential operators are 
the absolute values of the vector differences of the respec- 
tively indexed particles. Inherent to this ansatz the inter- 
nal degrees of freedom of the nuclei are not taken into 
account. But the screening forces of the electrons act at 
distances far beyond the range of the nuclear forces. The 
connection to the nucleus takes place via the calculated 
screening energy U e in the cross-section. Just as little any 
spin dependent effects are considered. So polarization ef- 
fects as observed in [19] cannot be described. Since the 
impact proceeds at keV-energies the inclusion of relativis- 
tic effects is not required. Furthermore, the velocity of 
the nuclei is small in comparison to the electron^. Hence, 
the Born-Oppenheimer approximation [20] can be applied 
which leads eventually to a drop out of the nuclear oper- 
ators in |[2|) , leaving ([21 H]) . The electron wave functions 
become then continuous functions of the nuclear coordi- 
nates. The movement of the projectiles is treated classi- 
cally owing to the acting forces. The now purely electronic 
Schrodinger equation can nevertheless only be solved ap- 
proximately, which is done based on the Ritz theorem by 
minimization of the energy functional E [<P] = (tf', bi<P) 
with \P from a restricted function space T' under control 
parameters. This is performed with the self-consistent field 
method by Hartree and Fock (HF). The electronic Hamil- 
ton operator is splitted up in a single electron operator 
((HI and a part containing the interaction potential (0]). 
The single electron problem has likewise solution functions 
ip (r, a) := ip (r) \ (c) with the space function ip and the 
spin function \- The solution of the multi-electron prob- 
lem is constructed with the anti-symmetrization operator 



2 Notations according to [18J. 

3 The electrons of the inner shells of heavy atoms attain rela- 



as a Slater determinant W = A <p a (1) • (fb (2) • ... • ip n (N) 
in accord with the Pauli principle. The variation of the en- 
ergy functional yields the Hartree-Fock equations [21] 
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with the Coulomb operator Jb and the exchange operator 
Kb defined as 

JbV'a (i) = [1pb (i) , — Ipb (j) J Ipa, (i) 



KbVa (<) = (j) , — ^ (j)) ^b (<) (6) 

They have the form of an eigenvalue equation for ip a with 
the orbital energy e a , though they are coupled by the op- 
erators Jb and Kb. Indeed, they are a system of coupled 
integro-differential equations which can only be solved it- 
eratively. To simplify matters the case of the restricted HF 
method is described where all n/2 orbitals are occupied 
with electrons. Otherwise the unrestricted HF method has 
to be used [22(123] . The HF method furnishes merely a 
stationary solution so far. A time development can be 
achieved with the unitary time propagation operator tfj (t) = 
U(t,t )ip(t ) which is given by U(t,t ) = e^*-* )" if 
the Hamilton operator is not explicitly time dependent 
[24]. Caley's form provides a unitary approximation for 
this equation as a calculable time discretization resulting 
in the time dependent Hartree Fock method using Hhf 
[25,26]. Such has been carried out by [27|I17] implement- 
ing the discretization of the electron wavefunctions on a 
space lattice. 

The HF method cannot include immediate electron- 
electron interaction^]- For electrons of equal spin the main 
correlation effects are incorporated in the Slater deter- 
minant. But the motion of electrons of opposite spin re- 
mains uncorrelated. To address this, several procedures 
have been developed which either start from the HF so- 
lution functions applying corrections or modify the HF 
method [23,28,29,30]. They are based on the configura- 
tion interaction method (CI) or the M0ller-Plesset per- 
turbation theory (MPn). The CI method uses the HF or- 
bital function and adjoins further wherein occupied or- 
bital functions are replaced by unoccupied virtual orbitals 
successively thus forming a series [22]. The linear coeffi- 
cients are determined by minimization of the energy func- 
tional. This series is infinite and would deliver an exact 
solution for the Schrodinger equation I13I4H in the limit. 
In practice, the series is aborted after the first to fourth 
term with few virtual orbitals in the individual terms. The 



4 The operators Jb and Kb in merely provide an averaged 
effective electron potential in the iteration due to the integra- 



tivistic velocities. This is commonly considered by corrections. tion in the scalar product, e.g. [18] 
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MPn method is essentially an implementation of the com- 
mon quantum mechanical perturbation theory [31] , It is 
usually aborted after the second to eighth order. Another 
approach is the density functional theory (DFT) [32"ll33]. 
Following the theorem of Hohenberg and Kohn a unique 
functional exists which provides an exact solution of (|3I4[| 
[34"ll35]. But there is no algorithm that can generate this 
functional so one needs to settle on heuristic ones and 
therefore the DFT is no ab-initio method in the end [36, 
[371l38ll39lHrjlliTlH2] . The DFT modifies the HF equations 
with an exchange interaction functional, which is to be 
supplied. 

The just described configuration interaction procedures 
can increase the computational costs by 1 to 3 orders of 
magnitude relative to the plain HF method O chapter 
6] . On the other hand the implementation of the orbital 
functions on space lattices is very cost-intensive and only 
feasible for few electrons. Hence the analytic HF method 
of Roothaan and Hall is used for larger systems [44,45, 
I46j . There the orbital functions are developed as a linear 
combination ij} a = J2 n C «A °f a se t of analytic func- 
tions {u n }l G L 2 (R 3 ) which are the basis of the finite 

dimensional product space T' = ® ^ T' n . Thus the HF 
equations are transformed to a system of pseudo-linear 
equations for the linear coefficients which are solved iter- 
atively again. The main computational costs now results 
from the integrations of the numerous scalar products. 
The computationally most efficient basis functions for it 
are the Cartesian Gaufi type orbital functions [17] 

u (GTO) = x i y j z k e-° tr2 , i,j,keN (7) 

with the sum I = i + j + k yielding the angular momen- 
tum quantum number of the orbital. Thereupon basis sets 
of increasing complexity are constructed being geared to 
actual atomic orbitals [23j. The computer resource usage 
scales with the total number of GTO basis functions as 
O (n) ~ n A customarily spanning two orders of magnitude 
|43l chapter 5]. The common minimal basis set is STO-3G 
where the atomic orbitals are approached by Slater type 
orbitals which contain the spherical harmonics and are 
constructed from 3 GTO's with fixed coefficients, there- 
fore called contracted [48|i49| . Split valence basis sets pro- 
vide two or more sizes of basis functions for the valence 
orbitals (e.g. D95 [50], LanL2DZ [51H52H53] . 3-21G [54ll55l 
[56], 6-31G [57 ll58ll59ll50ll6T] . 6-311G). Polarized basis sets 
allow the atomic orbitals to change their shape by pro- 
viding orbitals with angular momentum higher than for 
the ground state description, like p and d orbitals (here 6- 
31G(d,p) |62,63j). High angular momentum basis sets add 
more orbitals of p, d, and f type (here 6-311+G(3df,3pd) 
[64,65,66,67,68])- The inclusion of diffuse functions adds 
larger size versions of s and p type functions allowing the 
orbitals to cover larger regions of space which is advanta- 
geous for charge transfers in ionic bonds (indicated by a 
'+'). Most of the basis sets are only available for the lower 
periods of the periodic table. Only the LanL2DZ covers 
heavier atoms where the inner shells are described by ef- 
fective core potentials including relativistic effects. The 
extra basis (+XB) adds diffuse functions. 




Fig. 1. Cubic space centered crystal cell of LiH. 

In consequence, it becomes clear that the simulation 
of the pre-reaction impact inside the metal lattice involves 
so many electrons that the task has very high computer 
resource requirements. For the given workstation (Ultra- 
SPARC), the problem is actually oversized. So it is un- 
evitable to drop the time dependence and settle for the 
analytic HF method with restricted basis sets and few 
atoms, possibly augmented by the DFT. The impact is 
then treated as quasi static. The Coulomb adjusted rel- 
ative electronic force between the projectile and target 
nuclei is recorded and integrated along the trajectory T 
analogous to [27]117] yielding the molecular screening en- 
ergy 

U moi = J F re i • dr (8) 

T 

At the classical turning point the molecular screening en- 
ergy is identified with the screening energy U e = —U mo \ (r c i) 

So even for stationary orbital functions they are only 
feasible for relatively small molecules or require super com- 
puter scale power. 



3 Calculations 

The here adopted crystal structure data are taken from 
[69] . LiH has a cubic space centered structure with a = 
4.085 A [69l III, al] fig.[TJ Ta 2 H has a tetragonal structure 
with a = 3.38 A, b = 3.41 A 69 IV, j2] fig. H 

The calculations were performed with different crys- 
talline molecules, basis sets and models for the electron ex- 
change interaction employing the quantum chemical pack- 
age Gaussian [70]. One such molecule for LiH is depicted 
in fig. 02 The behaviour of the electrons in this molecule 
becomes clear by means of the density plots in the inter- 
section through the 4 Li atoms in fig. 01 The total propa- 
bility density of the electron wavefunction of the molecule 
is plotted in fig. [Ha. The maximum values are numeralized 
at the positions of the nuclei. The molecular wavefunction 
is not normalized to unity but to the number of electrons, 
i.e. L 3 <&*<]/ dr = n. So the probability density values can 
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Fig. 4. Different electronic densities in the 



be larger than 1 and may allow for certain statements 
on the electron number distribution. The electrons pref- 
erentially sojourn at the three fold charged Li nuclei, of 
course. Hence, the total molecular probability density is 
not very meaningful in order to observe the shift of the 
electrons due to the chemical bond. Much more suitable is 
the bond density (fig.[Hb) which is the difference between 
the total molecular probability density and the probabil- 
ity density of the single atoms. The result is a plot whose 
average value is zero and which shows the electron migra- 
tion when the bonds are formed. Now, it is visible that 
indeed an electron transfer towards the central hydrogen 
atom at the expense of the Li atoms took place. The den- 
sity matrix of the bond density is composed of two parts: 
the atomic net density (fig. [He) which results from set- 
ting the interatomic overlap equal to zero. One recognizes 
where the electrons are deducted from the atoms in order 




Li: -0.0523 H: 0.0343 
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Li: -0.0229 H: 0.0432 

molecule lih-lil with the basis 6-31G(d,p) 

to form the bonds; from Li in a stronger proportion than 
from H. The respective maximum values are marked above 
the plot and the minimum values below. In the overlap 
density (fig. [Hd) the atomic net density is set to zero and 
one recognizes the share of the molecular bond. Due to 
the bond formation the hydrogen atom profits by gaining 
electron shares. The ionic character of the bond becomes 
clear, since the overlap at the position of the hydrogen nu- 
cleus becomes a maximum and not in between the atoms 
as in the covalent bond. 

The projectile trajectory is chosen so that the projec- 
tile does not close in on other nuclei. According to [T7] on 
a nearby passage the spectating nuclei exert a repulsive 
force on the projectile causing a reduction of the screen- 
ing energy, on average. On the other hand ions traversing 
crystal lattices are deflected by the lattice atoms on tra- 
jectories inbetween them, known as the channeling effect. 
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Fig. 6. Relative force and molecular screening energy at the 
impact in the molecule lih-lil with the basis 6-31G(d,p). 




Fig. 2. Tetragonal crystal cell of Ta2H. 




Fig. 3. LiH molecule (lih-lil) with projectile trajectory and 
intersections. 



This means here that the crystal lattice of the metal atoms 
focuses the incoming deuterons on the interlaced lattice of 
the deuterons in the metal hydride [7.1 J counteracting the 
reduction of the screening energy in the D2 molecule of 

D21- 

During the quasi static simulation of the impact pro- 
cess instant recordings of the bond density at various dis- 
tances between the projectile and the target were taken 
which are assorted in fig. [5l The plots of the bond density 
lie in the intersection of the trajectory of the projectile 
which is drawn inclinated in fig. [3] and contains two Li 
atoms. For a contrast amplification the bond densities are 
not rendered with equal scale. At first the bond density is 
distributed at the projectile and the target. At a distance 



of 1.269 A an even enhanced transfer from the target to 
the projectile and the environment can be observed. From 
there on the electron transfer from the Li atoms to the hy- 
drogen nuclei approaching each other increases evermore. 
The calculations were performed for 100 points on the tra- 
jectory. The relative forces at these points are plotted in 
fig. [H The six points of the bond density plots (fig. \5§ 
are highlighted by rhombuses and mark interesting loca- 
tions. Beginning at a distance of 6 A the absolute value of 
the Coulomb-adjusted relative force raises up to a maxi- 
mum value at ~ 0.25 A and declines till zero for further 
decreasing distances which is understandable because the 
approaching nuclei virtually form a single one in compar- 
ison to the wavelength of the electrons. In so far the be- 
haviour is in accord with the TDHF calculations in [i~7[ fig. 
8] except for the oscillations due to the time dependence. 
The full curve in fig. [6] is the molecular screening energy 
which was calculated according to along the trajec- 
tory at every point. The final result is with —43.39 eV way 
too low compared to the experimentally obtained values 
(< 150 eV for L0 [H| and [190, 320] eV for the other met- 
als fTppp. 

The reasons for this unsatisfying result are to be sought 
in the hierarchy of model assumptions, approximations, 
and simplifications. Only very few of them could be reme- 
died within the scope of this work. However, some experi- 
ments with the bases could be employed whose restrictive 
impact on the movability of the electrons particulary in 
the conduction band is of jutting importance. The results 
of these calculations are summarized in the overview in 
fig. [7J They are arranged according to increasing complex- 
ity of the bases, in each case. Though from the smallest 
to the largest basis the value of l7 mo i becomes larger only 
by 10%. Yet the tendency is clear: the larger the basis the 
higher the screening potential. With an extension of any 
kind no physical error can be committed since the elec- 
trons arrange due to the energy minimization. Only the 



5 The results for the highly reactive alkaline metals Li and 
Na are heavily impaired by the surface oxidation under ion 
irradiation |11| . 
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Fig. 5. Instant recordings of the bond density during the impact in the molecule lih-lil with the basis 6-31G(d,p) 
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Fig. 7. Overview of the calculation results. 
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Fig. 8. Bond density in the molecule lih-lil with the basis 
STO-3G. 



restriction of the bases induces errors. Because the estab- 
lished available bases are centered at the atoms they are 
hardly suitable for the simplified description of the delo- 
calized electrons of the conduction band. 

The effect of the bases on the bond density of the 
molecule can be assessed by comparison with the smallest 
basis STO-3G whose bond density is depicted in fig. [H 
Unlike the 6-31G(d,p) basis in fig. [4jb no electron trans- 
fer to the hydrogen can take place here because the too 
small basis does not permit enough mobility for the elec- 
trons. Instead only between the nuclear positions a weak 
conglomeration accrues like for the covalent bond. 

The impact of the size of the molecule can only be 
estimated with narrow limitations, since by the addition 
of further spheres of H and Li atoms around the central 
target the problem would not be feasible for the worksta- 




Fig. 9. LiH molecule (lih-h2e) 



tion any more. Hence, only on one plane further H atoms 
were added as fig. [9] shows. The belonging density plots in 
the plane of the Li atoms (fig. LTO]) differ perceptibly from 
those of the simpler molecule (fig. [4j) . From the compari- 
son of the bond densities, the atomic net and the overlap 
densities it is apparent that the perimeter H atoms at- 
tract electrons on the expense of the central H atom and 
thus have a higher influx. Accordingly the screening po- 
tential is also lower (fig. [7J) . This behaviour is vitally due 
to the circumstance that the perimeter H atoms have no 
Li partners which supply their electrons. If those Li atoms 
were present the transfer of electrons from the central H 
atom would be quite perspicuously lower. It is interesting 
that the maximum values of the bond density form a ring 
around the central H atom and the maximum values of 
the perimeter H atoms are not centered around them but 
shifted towards the central H atom. 

Indeed the more complex basis brings higher values for 
the screening potential here again. Furthermore, it is here 
also documented how a method affects going beyond the 
Hartree-Fock method videlicet the density functional the- 
ory with the three parameter functional of Becke, B3LYP 
[38PTU39] . The more accurate provision for the electron 
correlation effects causes a stronger repulsion of electrons 
with equal spins and enhances the screening potential. 
However, this effect is not very large. 

A further possibility for a comparison is provided by 
the Mulliken population analysis which assigns the atoms 
of the compound a charge according to the electron trans- 
fer. Table [l] summarizes those results. The topmost atom 
is the central target H atom followed by the innermost 
sphere of metal atoms and further perimeter H atoms. 
The charges are specified in multiples of the elementary 
charge. The larger basis for Li provokes an almost ten 
times higher charge transfer to the H atom compared to 
the smallest one which yields almost 0.5 e and underlines 
the ionic character of the bond. The extended molecule 
with the four perimeter H atoms permits recognition of 
the reduction of the charge at the central H atom and the 
carry-over to the perimeter H atoms. In comparison the 
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Fig. 10. Different electronic densities in the 

rather high differences in the charge transfers were not 
reflected to this extent in the differences of the screening 
potentials (fig. [7|). The reason for it is that the charge 
assignment to the atoms is not a quantum mechanical ob- 
servable and hence tainted with a certain arbitrariness. 
This can also be seen at the smallest molecule whose Li 
atoms are fully symmetrically arranged where two of them 
donate and four sparsely ingest electrons. 

In the right part of the table charge assignments for 
Tantalum are listed. Owing to the considerable computer 
time consumption by virtue of the many electrons to be 
included only one simulation run on the smallest molecule 



the electron screening effect in metal deuteride crystals 
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molecule lih-h2e with the basis 6-31G(d,p) 

with the simpliest basis could be performed. The result 
of U e = 39.72 eV even stays behind the smallest Li value, 
however the charge assignment to the central H atom is 
not much smaller. The molecule and the projectile trajec- 
tory are shown in fig. [TTJ The behaviour of the relative 
force in fig. [32] is remarkable. The oscillations in the far 
field were caused by the large electron clouds with high 
angular momentum of the Ta atoms. 



For details refer to [3]. 
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Table 1. Atomic charges originating from the Mulliken popu- 
lation analysis. 
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Fig. 11. Ta2H molecule with projectile trajectory. 




Fig. 12. Relative force and molecular screening energy at the 
impact in the molecule Ta2H. 



4 Conclusion 

The calculated values for the electron screening energies 
stay clearly below the experimental ones. The simplifica- 
tions in the numerical model were enforced by the avail- 
able computing power. The performed numerical simula- 
tions consumed already a netto computation time of 1.6 a. 
Single points on the trajectory took one week. Anyhow, 
the simulations could show the migration of electrons from 
the host metal atoms to the hydrogen both for the target 
nuclei and the projectile in the impact. This migration 
becomes larger when the number of base functions is in- 
creased. However, in the electron transfer only the imme- 
diate adjacent atoms are involved and no distant atoms. 
This is due to the limited size of the basis which are de- 
signed for molecules and not for larger solids. The further 
extension of the basis sets would provide more freedom of 
mobility for the electrons. Particularly, the modelling of 
conduction band states as provided by Bloch wave func- 
tions could be an asset. But according to the computa- 
tional cost function already this extensions to the quasi- 
static model with an accompanying increase of the crystal 
size would require the next scale of computational power 
as computing clusters or super computers. Furthermore, 
time dependent effects can be expected to contribute the 
larger share in attaining the experimental screening results 
particularly when collective interactions between the col- 
lision partners and many electrons from a larger region 
are to be included like plasmons, convoy electrons and ef- 
fects of the Fermi surface. But a prerequisite for it are a 
high electron mobility provided by large basis sets. This 
would require the complete redesign of a computer code 
for a massive parallel computer worth many man years of 
work. But then very instructive insights into the electron 
dynamics for the understanding of the mechanism can be 
expected which is not accessible in the experiment so far. 
In exchange for higher code development effort the com- 
putational costs could possibly be reduced by three orders 
of magnitude using the new ^-operator formalism for the 
solution of partial differential equations by |72j . 
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